Mark A. O’Neill 


Faster Than Fast Fourier 


The fast Hartley transform is twice 
as fast as the fast Fourier and uses 
only half the computer resources 


J 


descendant. 

However, the fast Fourier transform (FFT) requires a large 
amount of computer resources. The fast Hartley transform 
(FHT) can accomplish the same results faster, using fewer re- 
sources. First described by Ronald Bracewell (see reference 1), 
the FHT handles many of the jobs now done by the Fourier 
transform. 

Both the FHT and the FFT let you map a continuous signal 
over time onto a frequency function. The Fourier transform 
maps a real function of time X(¢) to a complex function of fre- 
quency, F(f). 

The Hartley transform maps a real function of time X(t) onto 
a real function of frequency H(f). Since the Hartley frequency 
function is real, you need only single arithmetic operations to 
compute it. 

Compare this to the many arithmetic operations required of 
the complex Fourier frequency function (four operations for a 
complex multiply or divide, and two for complex addition or 
subtraction). 

Furthermore, real data arrays require only half the memory 
storage of complex data arrays. This means that the Hartley 
transform will require considerably less memory for a given 
data set than the Fourier transform. Therefore, the Hartley 
transform will be distinctly faster and use less memory than the 
conventional Fourier transform in digital filtering and image en- 
hancement applications where you have to process large 
amounts of data. 


oseph Fourier left to mathematics a rich legacy: 
The transform that bears his name processes 
many of today’s audio and electromagnetic sig- 
nals. Image processing and digital filtering of 
signals use the Fourier transform or its faster 





A Definition of the Hartley Transform 

Equation (1) shows the analytical form of the Hartley transform, 
and (2) shows its inverse transform, used to map the frequency 
function back into the time domain: 


A(f) = = \ x(t) cas(27ft) dt, () 


— oo 
oo 


X(t) = \ A(f) cas(2aft) df, (2) 


—-@ 


where cas(27ft) = cos(27ft) + sin(27ft). 


The function cas(27ft) was introduced by R.V.L Hartley (see 
reference 2), who first proposed the Hartley transform in 1942. 

You can see that these equations are very similar to those of 
the Fourier transform and its inverse: 


FN = a2 | xmeras 3) 


— oo 
oo 


\ F(f) e228 df (4) 


—-o 


where e/?7 = cos(2zft) +) sin(27ft), and e?"* = cos(2zft) 
— j sin(27ft). (These are known as Euler’s formulas. ) 

Note that I’m using the electrical engineering convention of 
labeling the imaginary unit i as j. The principal difference be- 
tween the two functions is that the real function cas(2 7ft) in the 
Hartley transform replaces the complex exponential term e*/?™* 
in the Fourier transform pair. 

The Hartley and Fourier transform functions in (1) through 
(4) deal only with continuous variables. As is often the case with 
computer data in the real world, signals are sampled at discrete 
intervals of time or for a specific interval. Fortunately, you can 
define a discrete transform that can represent a quantized con- 
tinuous signal, or a signal of limited duration. The discrete 
forms of the Hartley transform pairs are 


X(t) = 


N-1 
Pee 
H(f) = = ye Fo cas(2rft/N) (5) 
and ae 
N-1 
X() = J H(f) cas(2aft/N). (6) 
f=0 
Again, note the similarities to the Fourier transform pairs: 
N-1 
= al (-janfe/N) 
N=, »p Oe Sy, (7) 
t= 
N-1 
X() = DO F(f) etantr at, (8) 
f=0 
The Fast Algorithm 


As it stands in (5) and (6), computation of the discrete Hartley 
transform presents an analogous problem to the computation of 
the discrete Fourier transform. That is, you have to perform N? 
arithmetic operations to compute the discrete Hartley transform 
of an N-element data set. 

A classic paper by Cooley and Tukey (see reference 3) in 
1965 led to the development of a fast algorithm for the machine 
computation of a complex Fourier series. 

Essentially, the FFT uses a permutation process to bisect the 
data until data pairs are reached. Calculating the Fourier trans- 

continued 


APRIL 1988 * BYTE 293 


FASTER THAN FAST FOURIER 


form of such data pairs is trivial (1.e., rapid). 

The idea behind the permutation process is that it’s faster to 
split the data into pairs, compute the transform of the pairs, and 
recombine these to make the entire transform rather than to 
compute the transform for the complete data set. 

Permutation is particularly fast when the amount of data is 
large. If you superimpose all such two-element pairs using a 
process sometimes referred to as the “butterfly” (due to the ap- 
pearance of the diagram of the data flow; see figure 1), you can 
compute the Fourier transform of the input data set. It takes ap- 
proximately N log() seconds to compute the transform of an N- 
point data set. 

Bracewell (see reference 4) has shown that you can employ a 
similar methodology in the case of the Hartley transform. 
Again, you use the permutation process to bisect the data until 
you get data pairs. 

The Hartley transform of a data pair (a, b) is ‘’A(a + b, a — 
b), and the computation of such pairs is trivial. You can also 
superimpose these two-element sequences to calculate the Hart- 
ley transform of the input data set. However, to do so requires a 
formula that expresses a complete discrete Hartley transform 
(DHT) in terms of its half-length subsequences. 

Bracewell has shown by application ofthe Shift and the Simi- 
larity theorems that (9) expresses the general decomposition 
formula for the DHT. This general decomposition formula gen- 
erates the desired DHT by bisecting the data. 


Figure 1: The butterfly 
diagram for the fast Hartley 
transform for an eight- 
element data set. In the first 
column of cells, the indexes 
show the effects of 
permutation on the data 
elements. The next three 
columns represent the data 
elements as superimposed 
by the decomposition 
formula for the butterfly 
computation. These cells 
show the three indexes of the 
elements. You can see the 

ef fects of retrograde indexing 
for an eight-element data 
set in indexes for the third 
and fourth iteration of the 
butterfly. Radians at the 
bottom of the diagram give 
the angular spacing between 
trigonometric coef ficients. 
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Put another way, it’s the rule used to generate the elements to 
be used in the butterfly computation of the transform. You can 
apply similar methodology to the Fourier transform to yield the 
decomposition formula given in (10): 


H(f) = A,(f) + Hf) cos(2af/N.) 

+ H{N,—f) sin(27f/N,), (9) 
F(f) = Fi(f) + Fa( fye27t/™s, (10) 
where N, is the number of elements in the half-length sequence, 
and thus N,=N/2 for a data set of N elements. Figure 1 shows a 
complete butterfly diagram for computing the FHT for an eight- 
element data set. 

The decomposition formula for the FHT differs from the FFT 
in One important respect: The elements multiplied by the trigo- 
nometric terms are not symmetric. In the FFT decomposition 
formula of (10), terms multiplied by the trigonometric coeffi- 
cients involve terms only in F(f). In the FHT decomposition 
formula of (9), both H(f) and H(N,—f) have sine coefficients. 
This asymmetry becomes apparent when you express the dis- 
crete transforms as matrix operations: FFT matrix terms are 
symmetric about the matrices’ leading diagonal, while corre- 
sponding terms for the FHT are asymmetric. 

This introduces some computation problems, because asym- 
metric matrix processing is difficult to implement. You can deal 
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with this asymmetry by using an independent variable as an in- 
dex for the elements multiplied by the sine coefficients. This 
index decreases while the other indexes increase; this behavior 
is called retrograde indexing. Bracewell gives descriptions of 
the FHT and FFT algorithms using matrix formalism (see refer- 
ence 5). 

You can obtain the inverse Hartley transform by applying the 
FHT algorithm again to its own output, thus regenerating the 
input data. This means that you can use the same program code 
to compute the transform and its inverse. However, there is a 
slight asymmetry between the FHT and its inverse. In the case 
of the time-to-frequency transform, you need to scale the result 
of the butterfly computation. That is, for an input data set of V 
elements, you must divide the output of the butterfly by N to 
obtain the discrete Hartley transform. 

You don’t need to do this for frequency-to-time transform; 
consequently, the butterfly computation itself constitutes the 
inverse transform. It is not difficult to add a small amount of 
code to control whether scaling should be applied during a 
conversion. 


Comparing the FFT and FHT Algorithms 

You can use the FHT algorithm for many of the applications for 
which you would now use the FFT algorithm. These include 
convolution and deconvolution (used to remove artifacts intro- 
duced into data by imperfections in the sensors), and the genera- 
tion of power spectra for filtering. You can also obtain the Fou- 
rier transform itself from the Hartley transform. In fact, it is 
often faster to generate the Fourier transform and power spec- 
trum with the FHT than with the FFT, because computing the 
butterfly using real rather than complex quantities requires 
fewer floating-point operations. You assemble the real and 
imaginary parts of the FFT at the end of the calculation using 
the equations 


F, = A(f) + A(N-f) 
and 
Pin — A({) AW —)), 


where F;, is the real portion of the complex Fourier transform, 
F 18 the imaginary portion, and N is the number of elements in 
the data set. 

You can calculate the power spectrum directly from the Hart- 
ley transform using the equation 


Pf) = (HY + H(WN—-f)?)/2, 


where P, 1s the power spectrum. 

The theorem for convolving a pair of functions 1s almost iden- 
tical, whether you are considering the Hartley transform or the 
Fourier transform. Again, the FHT should prove superior to the 
FFT in terms of speed for any given implementation. Equation 
(11) summarizes the theorem for the convolution for the Hartley 
transform, (12) for the Fourier transforms. 


A@M® AO = Ai(f)A2-f) + Ai(—f)A2Af) 
AO® AM = Fi(f)F.f). 


The ® symbol denotes the convolution operation. 

The subscripts e and o in (11) denote the even and odd parts of 
the Hartley transform. Note that if one of the functions being 
convoluted is either even or odd, then the form of the convolu- 
tion theorem for the Hartley transform reduces to the particu- 
larly simple form indicated in (13): 


fi0® fp = A f)H.(f). 


(11) 
(12) 


and 


(13) 


In practical applications, many of the convolution functions 
are even. For example, the Gaussian function, used in image- 
enhancement work, is an even function. We can take advantage 
of the nature of these functions to use the computational shortcut 
provided by (13). 


Spectral Leakage 

As with the FFT algorithm, the FHT algorithm will produce an 

erroneous frequency function if the data set to be transformed 

does not smoothly approach 0 at both ends of its range. Such 

spectral leakage is undesirable in many cases. You can reduce it 
continued 


Listing 1: The code fragment from hartley.pas that 
shows the complete Hartley transform as implemented in 
TML Pascal. 


Fast Hartley transform routine 
version 1.00, Dated 22nd November 1986, 


transform:= forward; Forward transform from 
time frequency domain. transform := reverse; 
Reverse transform from frequency to time 
domain. 


power _ index: Index to which 2 must be raised 


to generate a transform containing ‘syze' 
elements: 
if syze = 8 then power index = 3; 
if syze = 16 then power index = 4; etc 
syze: Number of element in the input data 
array. 
Baan ca a Oa aaa a oe le a Se } 
Const 
datasize = 512; 
type 


direction_type = 
data_array type = 
real; 
var 
dir, test option: ehar; 
i,j,syze,iter,demo: integer; 
data array: data array type; 
Cransform direction: direction type; 


(forward, reverse) ; 
array[1l..datasize] of 


procedure fht(var data array: data_array type; 
power index, 
syze: integer; 
Cransform:directsom type) ; 


var 


Erg ind, 

beg inc, 

power, 

C. a, 

f a, 

i temp, 

section, 

S Start, 

5 6nd: integer; 

sne,csn: array[1l..datasize] of real; 
accu: array[1..2,1..datasize] of feal; 
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by multiplying the data set to be transformed by a suitable win- 
dow function before computing the transform. 
All these window functions cause the data to smoothly ap- 


P tats tine. Thi ti ah 
ermutation routine i aiataed proach 0 at the limits of its range. Stigall, Ziemer, and Hudec 


reorders the data before the butterly 


eer ama te reer eae | a ee have reviewed the effects of these window functions on the 
eee Le ee eee ee } power spectrum (see reference 6). 

If you’re interested in some experimentation, I’ve provided a 
function permute(index: integer): integer; number of these window functions below. The triangular win- 
var dow is a low-quality window function whose primary advantage 

i,j,s: integer; is its speed. An example of its use is to filter noise from telem- 
= _— etry signals in real time. 
ie million label a Triangular window: 
for i := 1 to power index do g ; 
Begin = Wn) = 2(n+0.5)/n, 
s := index div 2; . WN-n-1) = Wn), 
i ¢= J + 7 + dpdex = s — s; 
index := s; wheren = 0,1, ...i1,..., N/2. 
end; 
permute := j + 1; The Hanning and Hamming windows feature better reduction 
end; in spectral leakage, but at the expense of speed. I’ve used these 


window functions to preprocess satellite image data before using 
the Hartley transform to correct the images for defects in the 


Calculate the trigonometric functions : 
camera optics. 


required by the FHT and store values. 
For a N p@int transform, the trigno= 
metric functions will be calculated Hanning window: W(n) = 0.5(1 —cos(2a(n +0.5)/N)), 


at intervals of Nths of a turn 


Sa a ee et ee a } wheren = 0,1, ...1,...,N—1. 
procedure trig table(npts: integer); Hamming window: W(n) = 0.54 — 0.46(cos(27(n+0.5)/N)), 
Const 
pi = 3.14159265; wheren = 0,1, ...i1,...,N—1. 
var 
Le amit ; : : ; 
eras pee Am reall The Blackman window is for data sets where the Hamming/ 
begin Hanning windows would not be suitable. 
angle := 0; 
omega f= 2 * pi / npts; Blackman window: W(n) = 0.42 — 0.5cos(2a(n +0.5)/ 
EOr 10.4981 “ta mbes da (N—1)), 
begin 
sne(i] := sin (angle) ; wheren =0,1,...i,...,N-—l1. 
esn(i] := cos(angle); 
angle := angle + omega; 


Po You can find additional information about spectral window 
al functions in the papers by Harris (see reference 7) and Nutall 
(see reference 8). 


Calculate the address of the retrograde The FHT program 

index for the sine term for the dual The demonstrator program hartley. pas described below com- 
place algorithm, if it is required ... putes the Hartley transform and inverse Hartley transform of an 
“yea eam or we Hck pen eh es N-element data set. The program was originally written in 


Acornsoft ISO Pascal for the BBC Microcomputer, but few ma- 


function modify(power,s start,s end,index: 
ne “= aa wae: chine-dependent Pascal extensions are used. 


integer): integer; 


ae The routine has also been successfully portedto a Prime 9975 
gin ne 
if (s start = index) or (power < 3) then minicomputer, running Sheffield Pascal, and a Sun-3/160 
modify := index workstation running International Standard Organization level 0 
else Pascal. 
modify := s_start + s_end - index + 1]; As an example of just how easy it is to port the fast Hartley 
end, routine to other machines, I have included the code for hart- 
ee ee ee eae 0 te ee ley.pas as written in TML Pascal for the Macintosh. 


The hartley.pas program simply generates a suitable func- 
tion to be transformed and then calls the fht procedure twice. 


procedure butterfly (trig ind,i 1,i 2,i 3: Thus, the program calculates both the Hartley transform and its 

integer); 7 Sp 827 = inverse transform. The procedure fht is the major routine in the 

begin program. Listing 1 is a code fragment of hartley.pas that con- 
aceu [t a,1 1) 2 aecult a,i-1) + tains the fht procedure. 

accu[f_a,i_2] * csn[trig_ind] + The code within procedure fht complies with the ISO level 1 

accu[f_a,i_3] * sne[trig_ind]; standard for Pascal, except that the program uses the underscore 

continued character in variable names. Consequently, you can lift the pro- 

= continued 
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Puig: ind s= trig ind + syze div 2; 
aeeui(t a, 229i t= acew[f a,i 1) + 
acem PE a, i 2] * esn ftrag ind] + 
ace@in|f a,i 3] * smé[trigq iad); 
end; 


Main program for the fast Hartley 
transform. 


trig table(syze); 
for i t= 1 to syzéido weeul[f a,permute(i))) 


data_array[i]; 


Start of the Hartley butterfly 
transform ... 


power index do 


section := 1; 
tro’ inc := syze div (power + power) ; 


repeat 
trg ind #= 1; 
S stant, R= section * power + 1; 
s_end >= (section + 1) * power; 
for k := to power do 
begin 
butterfly (tg ind, ),] + power, 
modify(power,s start,s end, ] 
+ power)); 
tng ind := trg imd + trq ine; 
7 <= j + I? 
end; 
j = j + power; 
sectdon = section + 2% 
mntal 7 > syze-; 
power := power + power; 
1 temp := t_a; 
Gu = £ a? 


End of Hartley butterfly. The results are 
scaled if necessary, and then placed in 
back into the array data ... 


case transform direction of 
forward: for i := 1 to syze do 
data_array[i] 
aceulf a,i)] / syae; 
reverse: for i := 1 to syze do 
data_array[i] 
aeéu (Ea, i); 
end; 
end; 
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cedure from the demonstration program and use it freely in your 
own applications. | 

However, a brief description of fht and its embedded subpro- 
cedures may be useful if you’d like to experiment with the 
program. 

The procedure trig_table precalculates the trigonometric 
functions required by fht. For an input data set containing N 
elements, adjacent entries in this table are 27/N radians apart. 
This precalculation of the trigonometric terms avoids redundant 
computations of these values if they are needed frequently dur- 
ing processing. 

The procedure permute bisects the input data set progres- 
sively until data pairs are reached. The algorithm given is 
adapted from one described by Bracewell, who has also de- 
scribed faster permutation routines (see reference 5) that you 
can substitute if you want. 

The procedure butterfly calculates the butterfly of a single 
index pair. This procedure requires three indexes to use as entry 
parameters. This takes into account the possibility of retrograde 
indexing for the data element multiplied by the sine factor. 
Where retrograde indexing is necessary, the procedure modify 
calculates the retrograde indexes used by the procedure 
butterfly. 

The main program of the fht procedure first precalculates 
the trigonometric function tables and permutes the input data. 
The program then enters an iterative loop that computes the but- 
terfly for the whole of the N-element input data set. You then 
scale the resulting output by a factor of N if you are computing 
the time-to-frequency transform. 

When you have to deal with the complex form of the fre- 
quency transform, you can obtain Fourier transform from the 
FHT algorithm by operating on the output data of the procedure 
fht using the Pascal procedure get_FFT.pas, shown in listing 
2. By the time you exit from this procedure, the real and imagi- 
nary parts of the Fourier transform will be in the arrays. They 
are passed via the dummy conformant array variables r_pt and 
i_pt, respectively. 

A similar procedure in get_pwr.pas called get_power_ 
spectrum (shown in listing 3) allows you to obtain the power 
spectrum from the Hartley transform. At exit from this proce- 
dure, the power spectrum will be in the array passed to the pro- 
cedure via the dummy conformantarray variable p_sp. Youcan 
use both of the above procedures in applications programs 
without any modification. These procedures conform to the ISO 
level 1 standard of the Pascal language, and should be usable 
directly on any computer that supports an ISO level 1 Pascal 
compiler. 

The modifications required to compile these procedures with 
an ISO level 0 compiler are trivial. Basically, it is just a case of 
replacing all the conformant array parameters ’ procedure head- 
ers of get_fft and get_power_spectrum, with array param- 
eters of fixed size. 


Performance 
The BBC Microcomputer implementation of the FHT program 
has proved to be fast and is currently running on a 6502A-based 
BBC Microcomputer with a second add-on 3-MHz 6502 pro- 
cessor, which executes the program. This hardware can com- 
pute a 256-point transform in about 32 seconds. This compares 
favorably with the 2 minutes it took Bracewell’s original imple- 
mentation of the FHT algorithm to accomplish a 256-point 
transform on an HP-85 system. You can expect even better per- 
formance from 16-bit microcomputers. 

Table 1 shows the timings for running the FHT program on 
the Macintosh Plus and Mac II using TML Pascal and the Sun- 
3/160 workstation using ISO Pascal. 


continued 
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Table 1: Timings for the hartley.pas demonstrator program on several machines. For the Macintosh computers, TML 
Pascal version 2.50 was used, and times are to the nearest second. No 68881 coprocessor code was generated for the test, 
so the Mac II’s 68881 was not used. You can expect better times if you use the 68881. For comparison, the results of the 
same code compiled on a Sun-3/160 workstation using ISO Pascal are shown. Times are in seconds. 


Mac ill 
FPU not used 


Sun-3/160 
FPU not used 


Sun-3/160 
68881 FPU used 


Mac Plus 
FPU not used 


Number 
of points 


2 0.00 
4 0.00 
8 0.00 
16 1.00 
32 1.00 
64 2.00 
6.00 
12.00 
27.00 


0.00 
0.00 
0.00 
0.00 
1.00 
1.00 
2.00 
3.00 
6.00 


0.00 
0.02 
0.04 
0.08 
0.18 
0.34 
0.70 
1.48 
3.60 


0.00 
0.00 
0.00 
0.02 
0.04 
0.06 
0.08 
0.18 
0.40 


Listing 2: Pascal code to compute the Fourier 
transform from the FHT algorithm. 


Precedure get FFI (var data: 
array(l0..hOss_int] of rem; 

Var r pe+ array(11..hDas got) of real; 
Var 1. pt) arrmayll2..n2es ime] of 1 éal; 
Size? & amici; 

Vam 22 Sime; 

begin 

1 = As 
while i<size do 

begin 

be PEli)  -= datalij + dataleize — 1. + 1]; 
1 peli] t= datai) = “datalisize = 1 + 1)¢ 
1 t= 1&4 2; 

end; 
end; 


Note that you can attain very high speeds with the Sun work- 
Station, with its combination of a 20-MHz 68020 processor, 
68881 floating-point coprocessor, and optimizing Pascal 
compiler. 


New Limits to Explore 

I hope this article has given you enough information to seriously 
consider using the Hartley transform for your signal processing 
needs. Since the transform uses only real functions, you don’t 
need computationally expensive complex math to digitally filter 
or enhance a signal. 

The elimination of complex numbers also reduces the amount 
of memory required to process a signal. Finally, since the Hart- 
ley transform uses fewer operations to process a signal, you will 
have fewer roundoff errors. 

The FFT made a lot of what we call image processing possible 
by manipulating large amounts of data in a reasonable amount of 
time. The FHT offers better performance using less computational 
resources. With the same code you use to compute the trans- 
form, you can compute its inverse when necessary. It will be 
interesting to see what new uses will result from the expanded 
limits of processing that the Hartley transform provides. @ 


Editor’s note: Source code listings that accompany this arti- 
cle, hartley.pas, get_FFT.pas, and getpwr.pas, are avail- 
able on BIX, on BYTEnet, on disk, and in the Quarterly Listings 
Supplement. See page 3. 
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Listing 3: Pascal code to compute the power spectrum 
from the FHT algorithm. 


procedure get _power spectrum 

(var datas array(1@..h@:s_ int] of real; 
Var B sp: areayirll..hlss.gnt| Of rea; 
Siz6: mS imt) > 

Car it Ss tae, 


while i<size do 
begin 
P splij) <= (date[i] * datali} + 
data[size - i+1] * 
data({size - iJ) / 2; 
1 f=] 2 tools 
end; 
end; 





REFERENCES 

1. Bracewell, Ronald N. ‘‘The Fast Hartley Transform.” Proceed- 
ings of the IEEE, vol. 72, no.8, p.1010, 1984. 

2. Hartley, R. V. L. “A More Symmetrical Fourier Analysis Ap- 
plied to Transmission Problems.” Proceedings of the IRE, vol. 30, 
p. 144, 1942. 

3. Cooley, J. W., and J. W. Tukey. “An Algorithm for the Ma- 
chine Computation of Complex Fourier Series.” Mathematical 
Computing, vol. 19, p. 297, 1965. 

4. Bracewell, Ronald N. “The Discrete Hartley Transform.” Jour- 
nal of the Optical Society of America, vol. 73, p. 1832, 1983. 

5. Bracewell, Ronald N. The Fast Hartley Transform. New York: 
Oxford University Press, 1986. 

6. Stigall, P., R. E. Ziemer, and L. Hudec. “‘A Performance Study 
of 16-bit Microcomputer-Implemented FFT Algorithms.” JEEE 
Micro, p. 61, November 1982. 

7. Harris, F. J. “On The Use of Windows for Harmonic Analysis 
with Discrete Fourier Transforms.” Proceedings of the IEEE, vol. 
66, no. 1, p. 51, 1978. 

8. Nutall, H. H. “Some Windows with Very Good Sidelobe Behav- 
iour.” IEEE transcripts of Acoustics Speech and Signal Processing, 
vol. ASSP-29, no. 1, p. 84, 1980. 


Mark A. O'Neill is a member of a research team working on the 
automated production of topographical maps from satellite 
(SPOT) images. He is also a member of the department of photo- 
grammetry and surveying at University College in London, U.K. 


